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Abstract 

We use molecular dynamics simulations of the SPC-E model of 
liquid water to derive probability distributions for water density fluc- 
tuations in probe volumes of different shapes and sizes, both in the 
bulk as well as near hydrophobic and hydrophilic surfaces. To obtain 
our results, we introduce a biased sampling of coarse-grained densities, 
which in turn biases the actual solvent density. The technique is easily 
combined with molecular dynamics integration algorithms. Our princi- 
pal result is that the probability for density fluctuations of water near a 
hydrophobic surface, with or without surface-water attractions, is akin 
to density fluctuations at the water-vapor interface. Specifically, the 
probability of density depletion near the surface is significantly larger 
than that in bulk. In contrast, we find that the statistics of water den- 
sity fiuctuations near a model hydrophilic surface are similar to that 
in the bulk. 



1 Introduction 

According to theory, fluctuations of water density on sub-nanometer length 
scales obey Gaussian statistics [H HI |3] , while those on larger length scales 
deviate significantly from this behavior [4^. The deviations at large length- 
scales reflect the fact that water at standard conditions lies close to water- 
vapor coexistence jH E]- A large enough repellent surface can therefore 
induce the formation of a water-vapor-like interface, and as such, the prob- 
ability of water depletion is enhanced near such a surface. In particular, the 
presence of this liquid- vapor-like interface facilitates large fluctuations in its 
vicinity compared to that in the bulk. This perspective is at the heart of 



current ideas about hydrophobic effects [6l [7] . In this paper, we use molec- 
ular simulations to further examine the validity of this perspective. We 
demonstrate marked similarities between water-vapor interfaces and water- 
oil interfaces for large enough oily surfaces with low radii of curvature. 

The fact that a purely repulsive hydrophobic surface induces the forma- 
tion of a vapor-liquid interface is well-established jH [3 [9l [10] . Dispersive 
attractions between the hydrophobic surface and water can mask this effect, 
as the attractions move the interface to a mean position immediately adja- 
cent to the hydrophobic surface. This removes any significant presence of 
vapor on average [TOl [TTl [121 [13] . As a result, it is difficult to detect the pres- 
ence of a water- vapor-like interface by considering only the mean behavior of 
water density profiles. Rather, the presence of this interface is more directly 
reflected in fluctuations away from those profiles [HI [H] [THl [13 [13 [IH] , 
which motivates our focus here on the statistics of these density fluctuations. 

Large length-scale water density fluctuations control pathways for assem- 
bling hydrophobically stabihzed structures [3l|2nil2lll22l|23l|2l[25l|26ll271. 
Several groups have studied the statistics of these fluctuations [3l [28] [29] . 
This paper is distinguished from this earlier work by the fact that we report 
the statistics of large length-scale fluctuations, which are expected to be 
fundamentally different from the small length-scale fluctuations studied in 
the earlier work. 

The central quantity to be examined is the probability for finding N 
water molecules in a sub- volume v, P^{N). Hummer and Pratt and their 
co-workers introduced the idea of studying this function as a route to under- 
standing solvation j3l |28l |30l |31]. In pure solvent, large length-scale fluctu- 
ations are very rare. These are the fluctuations of interest here. Therefore, 
to obtain reasonable statistical information, we must employ some form of 
biased non-Boltzmann sampling |32| |33] . The simulation setup and system 
details are described in the next section, after which we explain the specific 
sampling method that we employ. Results are then presented, beginning 
with those for the length-scale dependence of fluctuations in bulk water, 
followed by a discussion of the effect of dispersive attractions between the 
water molecules and a hydrophobic surface. Finally, we contrast fluctuations 
near a hydrophobic surface with those near a hydrophilic surface. 

Our results are consistent both with the Lum-Chandler- Weeks (LOW) [3] 
theory for the length-scale dependence of hydrophobic effects and with re- 
cent computer simulation studies |18| IM] . In this sense, our results may 
seem unsurprising. This paper, however, is the first publication of com- 
puter simulation results demonstrating expectations of the LOW theory for 
fluctuations that deviate far from the mean in an atomistic model of water. 
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2 Simulation Models 



Density fluctuations of SPC-E water [35j in the canonical ensemble at 298K 
can be quantified using the LAMMPS molecular dynamics (MD) simulation 
package. [36] Specifically, we evaluate the probability, P^{N), of finding N 
water oxygens in a probe volume of interest, v. A straight-forward canonical 
ensemble simulation, with an average water density, p = 1 g/cm^, would 
suppress large density fluctuations. To avoid this suppression, a particle- 
excluding field is introduced on one surface of the simulation box, with the 
box large enough that the net density is less than the bulk liquid density. 
Through this construction |37| [38] we ensure that the bulk liquid remains in 
the center of the box, that it is at coexistence with its vapor phase, and that 
a free liquid-vapor interface acts as a buffer in the event of a large density 
fluctuation. Solutes are placed near the center of the box, deep within the 
bulk liquid and far from the free water-vapor interface. 

To model a hydrophilic solute, we consider water molecules within a 
particular sub- volume of bulk water, specifically a sub- volume of dimensions 
3 X 24 X 24 A'^. Taking a configuration of equilibrated bulk liquid, we 
immobilize those water molecules, which are within that considered sub- 
volume. The immobilized molecules are then the solute. Since the surface 
of this idealized hydrophilic solute is made up of fixed water molecules in a 
configuration that is typical of bulk water, the surface is unlikely to disrupt 
the hydrogen bond network of the neighboring solvent. 

To model a hydrophobic solute, we replace the fixed waters in our model 
of a hydrophilic solute with methane-like oily particles. These oily particles 
are uncharged and interact with the surrounding solvent water molecules via 
a standard water-methane Lennard-Jones potential. The resulting potential 
energy field expels water from the volume occupied by the solute. It also 
attracts water with dispersive interactions. The solute constructed in this 
way and a typical configuration of solvent water are shown in Fig. 1. 

In order to study the effect of dispersive attractions on water-density 
fluctuations, the Lennard-Jones (LJ) pair potential between the hydrophobic 
solute and the solvent is split into repulsive and attractive parts using the 
Weeks-Chandler- Andersen (WCA) prescription. [5^ The role of attractions 
can then be examined systematically with a scaling parameter A, 

uxir) = uo{r) + XAu{r), (1) 

where uo{r) and An(r) are the WCA repulsive and attractive branches of 
the LJ potential, respectively, a = 3.905A and e = 0.118Kcal/mol were 
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used as the LJ parameters for the oily particles [IQ], and Lorentz-Berthelot 
mixing rules were used to obtain water-solute interaction parameters. 




Figure 1: A typical probe volume, v = (3 x 24 x 24) A^, is shown here 
adjacent to a hydrophobic surface (blue particles). The rendered water 
molecules (red and white) are in a typical equilibrium configuration taken 
from one of our simulations. 



3 Umbrella sampling 

Since we have chosen molecular dynamics to probe our system, it is conve- 
nient to use a biasing umbrella potential that produces continuous forces. 
The potential we use is a function of the entire set of water oxygen co- 
ordinates: {fj}, i = 1,2,...M. The number of molecules, A^, in a specific 
volume, V, is not a continuous function of {rj}, but the biasing potential we 
use to influence this number is a continuous function of these variables. In 
particular, we focus on the coarse-grained particle number: 

„ M 

iV({ri},v) = / dr^$(r-r,)/iv(r), (2) 

i=l 
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where h^{r) = 1 or 0, depending on whether or not r is in volume v, r has 
Cartesian coordinates x, y, z, and 

<I>(r) = 4>{x)<P{y)H^) (3) 

with (j){x) being a normahzed, truncated and shifted Gaussian-hke distribu- 
tion, 

^{x) oc [eM-x^/2f) - eM-rcf^fMrc - \x\). (4) 

The proportionahty constant is the normahzation constant, d{x) is the 
heavy-side step function, and we have chosen the coarse-graining length 
^ to be O.lA and the cut-off length Tc to be 0.2A. 

In the limit ^ ^ 0, the dynamical variable iV({rj},v) is the actual num- 
ber of water molecules in the volume v. But for finite ^, A^({rj},v) is a 
continuous and differentiable function of {rj}. We construct the biasing 
potential with this variable. In particular, we let 

C/({ra;K,r?) = ^[iV({r,},v)-r/]', (5) 

where k is a positive constant. We have found it convenient to use k = 0.25 
kcal/mol. Simulating our system in the presence of this umbrella potential 
allows us to bias the system towards configurations with A^({rj}, v) values 
near rj. When 7] differs substantially from the mean value for the dynam- 
ical variable, these configurations will be very improbable in the unper- 
turbed system. Nevertheless, these configurations can be accessed reversibly 
through a series of simulations that slowly change the control variable ry. 

By influencing the coarse-grained number of particles in the probe vol- 
ume, the control variable r] also influences the actual number of water 
molecules in that volume. We have picked a small value of the coarse grain- 
ing length ^ to ensure that the latter influence is significant. As a result, 
with a series of simulations with different values of r], histograms for both the 
coarse grained number and the actual number can be collected and then un- 
biased and stitched together within the framework of the weighted histogram 
analysis method (WHAM) |4H H2| SSj. This procedure yields the joint 
distribution function that the coarse grained particle number, iV({rj},v), 
has value N and the actual particle number is N. The joint distribution, 
Pv{N, N), can then be integrated to give the distribution of interest, P^^N). 
With P^{N) known, the free energy of solvation of a "hard" probe cavity of 
volume V, A/Uy, is also known because ^ 

/3A/Uv = -lnPv(O) (6) 

where /cb/? = ^/T is inverse temperature and /cb is Boltzmann's constant. 
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Figure 2: (a) Probability distribution of finding water oxygens in a probe 
volume V in bulk water for a "small" cubic v = (6x6x6) A^, a larger "cubic" 
V = (12x12x12) A3 and a "thin" v = (3x24x24) A^. The solid line refers to 
the Gaussian distribution with the same mean and variance; 6N = N—{N)^ 
refers to the instantaneous fluctuation in N from its mean, {N)^. (b) The 
solvation free energy, Afi^, in units of k^T, per unit surface area. A, for 
probe cavities with different thicknesses and square cross-sections, but the 
same large volume [(12A3) = 1728A3]. The dashed line indicates the value 
of the the surface tension of SPC-E water. 
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4 Results 



4.1 Length-scale dependence of density fluctuations in bulk 
water 

In Fig. 2, we show results for P^{N) in different probe volumes in bulk 
water. Because Pv{N) is Gaussian for molecularly sized probe volumes jH], 
we present these results in comparison with Gaussians of the same mean, 
(A^)v, and the same variance, {{5N)'^)^. For small deviations from the mean, 
P^{N) is essentially Gaussian, and for small volumes v, only small deviations 
from the mean are possible. For large v, however, in the wings of the 
distribution for small N, the distribution differs markedly from Gaussian. 
In pure water, the chance of observing these deviations is negligible, less 
than one part in many powers of ten. On the other hand, these deviations 
become accessible and even dominant near a sufficiently large and repellent 
solute particle. This is true because in these wings of the distribution the 
free energy to reduce N, namely —k-BTln [P^{N)], can vary linearly or sub- 
linearly with N. See Fig. 2a for the case of a large cubic probe volume. As 
a result of this behavior, the introduction of a potential energy, which scales 
linearly with the number of particles it couples to, can favor low N to the 
point where low values become the most probable values. It is this shift in 
the distribution, which can occur only for large v, that is responsible for the 
large length-scale hydrophobic effects 

In the case of a large but thin rectangular volume, the wings of P^^N) 
also exhibit a deviation from Gaussian behavior. However, in this case, the 
distribution lies below the Gaussian. Despite the fact that both the large 
cubic volume and the thin rectangular volume have the same volumetric 
size, the probabilities for emptying the cubic and the thin probe volumes 
differ by ~ 25 orders of magnitude. This behavior can be rationalized within 
the context of large length-scale solvation behavior being governed by the 
energetics of interface formation. The thin volume has a much greater sur- 
face area than the cubic volume. Indeed, the A/Uv computed from Eq. (6) 
for several probe volumes with fixed volumetric size but varying thickness 
confirms this rationalization, as shown in Fig. 2b. 

4.2 Water near hydrophobic surfaces and the effect of dis- 
persive attractions 

Figure 3 shows normalized mean densities as a function of the distance, 
X, from center of the idealized large flat hydrophobic solute (see Fig. 1). 



7 



1.5 




3 6 9 12 



x(A) 

Figure 3: Mean water density p{x), relative to its bulk liquid value, pb) per- 
pendicular to extended hydrophobic solutes with different levels of dispersive 
attractions quantified by A. 



Several strengths of dispersive attractions between the solute and water are 
considered. See Eq. (1). For A = 0, the mean density profile is sigmoidal, 
suggestive of a vapor-liquid interface. However, addition of a small amount 
of attraction results in a qualitatively different density profile. For A = 
0.4, there is a maximum in the density profile accompanied by layering. 
Further increasing the attractions leads to a more pronounced maximum and 
layering. This behavior is in accord with the predictions of LCW theory [10 . 
Nevertheless, contrary to LCW theory, it has been suggested that a layered 
density profile implies an absence of a liquid-vapor-like interface near an 
extended hydrophobic surface with dispersive attractions to water. [TT| [Hj. 
The source of this incorrect impression is that the mean density by itself is 
not an indicator of liquid-vapor-like interfaces. Interfaces are relatively soft 
so that a weak perturbation can affect the location of the interface and thus 
the mean density profile while not destroying the interface. In other words, 
in order to fully appreciate the effect that a hydrophobic solute has on the 
surrounding solvent, one has to look at not only at the mean density but 
also at the density fluctuations. 

The statistics for these fluctuations can be obtained from the distribution 
of particle numbers in suitably chosen probe volumes. Figure 4a shows 
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Figure 4: (a) P^{N) for probe volumes near hydrophobic solutes with dif- 
ferent A values, (b) The corresponding solvation free energies. For A < 3, it 
is easier to solvate v near the hydrophobic solute than it is to solvate v in 
bulk water. 
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Pv{N) distributions for the thin rectangular probe volume v=(3 x 24 x 24) 
placed between x = bA and x = SA. With this position, there is no overlap 
between the van dcr Waal's radii of solute particles and water molecules in 
V. The distributions for A = and A = 0.4 are similar, with the probability 
of density depletion slightly lower for the latter case, but still significantly 
higher than that in the bulk. 

From Fig. 4b, we see that the free energy to empty this probe volume 
adjacent to the large hydrophobic solute with A = 0.4 is blksT, whereas 
the free energy to empty this same v when it is in bulk and far from the 
hydrophobic surface is 147A;bT'. We can understand this large difference in 
terms of the free energetics of interface formation. As seen in Fig. 3, an 
extended hydrophobic solute with A = leads to the formation of an inter- 
face due to the unbalancing of attractions in water. Py{0) for v adjacent to 
the solute is then essentially the probability to move the interface outwards 
by SA from a;~5Atox~8A. The free energetic penalty, A/iv for this 
process is ASksT as shown in Fig. 4b. The corresponding A/iy for A = 0.4 
is 57kBT, only a 9kBT increase on turning the attractions on to A = 0.4 
and as much as QOk^T less than that required to form interfaces. Hence, 
while the presence of a density maximum and layering at A = 0.4 might lead 
one to question the presence of an interface, the fluctuations in density or 
the ease with which a volume near the hydrophobic surface can be vacated 
leave no doubt as to its presence. Furthermore, Fig. 4b shows that one has 
to ramp up the attractions to nearly 3 times the value of typical dispersive 
attractions for the solvation free energy near a hydrophobic solute to equal 
that in the bulk. 

4.3 Fluctuations near hydrophobic and hydrophiUc surfaces 

PviN) for the probe volume v = (3 x 24 x 24) A^, near the hydrophobic solute 
(A = 1) is compared with that for the probe volume near the hydrophilic 
solute in Fig. 5a. P-v{N) near the hydrophilic solute (circles) is nearly 
identical to the Pv{N) for v in the bulk, (squares) On the other hand, the 
probability of density depletion for v near the hydrophobic solute (diamonds) 
is significantly higher. In Fig. 5b, we show the solvation free energy, A/iv, 
of the probe cavity as it is moved away from the solute, as a function of 
the distance between the center of the solute and the center of the cavity 
volume. While A/Xv for v near the hydrophilic solute (circles) remains equal 
to that for v in the bulk (solid line), A/Xy for v near the hydrophobic solute 
(diamonds) increases monotonically as v is moved away from the solute 
and eventually plateaus to the bulk A/Xy- Our results indicate that the 
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Figure 5: (a) P^{N) for v = (3 x 24 x 24) near a hydrophobic solute, near 
a hydrophilic solute, and in bulk, (b) The change in solvation free energy of 
the probe volume v as it is moved away from the solute. 
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hydrophobic surface affects density fluctuations in the water at a distance of 
up to ~ lOA. These results are in agreement with recent simulation studies, 
reporting the free energy of solvating a molecularly sized WCA cavity near 
hydrophobic surfaces |il8j and also the potential of mean force for bringing 
two hydrophobic plates close together. [3l] 

These results bear directly on nano-scale assembly. In particular, we have 
shown that the probability of water density depletion near a hydrophobic 
surface is significantly enhanced. When two such hydrophobic surfaces ap- 
proach each other, at a particular separation, the liquid between them is 
sufficiently destabilized to make drying and hydrophobic assembly kineti- 
cally accessible. In contrast, density fluctuations near a hydrophilic surface 
are identical to those in the bulk and the vapor phase is not stabilized by 
the presence of a hydrophilic surface. Hydrophobic and hydrophilic surfaces 
thus differ fundamentally in the way they affect the fluctuations of water 
molecules in their proximity. It is not the mean density, but rather the 
statistics of fluctuations that is most important. 
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